create y 1973 2008
'importing data from Excel for Pakistan
import  "C:\Users\Maryam\Desktop\BS Studies\PhD Thesis-II\EViews and STATA Progarm Codes\Chapter-9\Chapter 9.xlsx" range="Pakistan"

'***************************************************************************************
'ESTIMATING BALASSA-SAMUELSON EFFECT FOR rer_def_NT, rop, gdppc, gexp & a_tilde
'**************************************************************************************
'*************************************************************
'STEP 0: Tests for Unit Root in Individual Time Series
'*************************************************************

'***************************************************************
'Graph for Pakistan's  rop
'***************************************************************
                                        
genr rop = rop
freeze(figure_rop) rop.line
figure_rop.addtext(t) rop (Pakistan):  1973-2008
figure_rop.addtext(b) Year
figure_rop.addtext(l) rop
figure_rop.legend(off)
                                                 
'We see from the FIGURE that rop has a time trend to it.  So we would include both an intercept and a time trend in our unit root regression equtions. 

'*********************************************************************************
'ADF Unit Root Test for Pakistan's rop
'We now run our first ADF test
'**********************************************************************************
 
freeze(table_9_8_rop_adf) rop.uroot(adf,trend,info=sic)

'Note that the SIC automatic lag selection picks lags, p = 0.  The unit root test produces a t-value of -2.23 which is greater than our 5% criterion -3.54.  Thus, at this point, we cannot reject the null of a unit root.

'Now, let's check for white noise.  To do that, I first set all the residuals = 0, then run the adf test and finally will check for white noise.

genr resid = 0
freeze(mode=overwrite,rop_adf) rop.uroot(adf,const,trend,info=sic)
freeze(rop_adf_correl) resid.correl
 
'Based on the Q-statistic, I conclude that the residuals are white noise.  Putting it all together, I conclude that the rop series is not level stationary.

'The next thing I do is test whether the differenced series is stationary using the ADF test.  I once again begin by graphing the (differenced) series.
 
genr ropdiff = d(rop)
freeze(figure_ropdiff) ropdiff.line
figure_ropdiff.addtext(t) Drop (Pakistan):  1973-2008
figure_ropdiff.addtext(b) Year
figure_ropdiff.addtext(l) Drop
figure_ropdiff.legend(off)

'The graph is not particularly illuminating.  Depending on how you look at it, it could have a time trend to it.  However, since the previous ADF test concluded that the equation for rop do have a time trend, we run the ADF test with both a constant as well as a time trend.

'So we begin the whole process over again: 

genr ropdiff = d(rop)
freeze(table_9_8_ropdiff1_adf) ropdiff.uroot(adf,const,info=sic)

'Note that the SIC roptomatic lag selection picks no lags, p =0.  The unit root test produces a t-value of -7.59 which is now smaller than our 5% criterion -2.95.  Thus, we may now reject the null of nonstationarity in first differenced series of rop.  There is no reason to go further.  The last thing we do is check for white noise.

genr resid = 0
freeze(mode=overwrite,ropdiff1_adf) ropdiff.uroot(adf,const,info=sic)
freeze(ropdiff1_adf_correl) resid.correl

''Based on the Q-statistic, I conclude that the residuals are white noise.  Putting it all together, I conclude that the rop series is I(1).

'*********************************************************************************
'DF-GLS Unit Root Test for Pakistan's rop
'We now run our first dfgls test
'**********************************************************************************
 
freeze(table_9_8_rop_dfgls) rop.uroot(dfgls,trend,info=sic)

'Note that the SIC automatic lag selection picks lags, p = 0.  The unit root test produces a t-value of -1.87 which is smaller than our 5% criterion -3.19.  Thus, at this point, we may reject the null of a unit root.

'Now let's see if the series is difference stationary or not

genr ropdiff = d(rop)
freeze(table_9_8_ropdiff1_dfgls) ropdiff.uroot(dfgls,const,info=sic)

'Note that the SIC automatic lag selection picks no lags, p = 0.  The unit root test produces a t-value of -3.94 which is now smaller than our 5% criterion -1.95.  Thus, we may not reject the null of nonstationarity in first differenced series of rop.  
 
''Putting it all together, I conclude that the rop series is I(1), a finding compatable with my ADF test results.

'***************************************************************
'Graph for Pakistan's  gexp
'***************************************************************
                                        
genr gexp = gexp
freeze(figure_gexp) gexp.line
figure_gexp.addtext(t) gexp (Pakistan):  1973-2008
figure_gexp.addtext(b) Year
figure_gexp.addtext(l) gexp
figure_gexp.legend(off)
                                                 
'We see from the FIGURE that gexp has a time trend to it.  So we would include both an intercept and a time trend in our unit root regression equtions. 

'*********************************************************************************
'ADF Unit Root Test for Pakistan's gexp
'We now run our first ADF test
'**********************************************************************************
 
freeze(table_9_8_gexp_adf) gexp.uroot(adf,trend,info=sic)

'Note that the SIC automatic lag selection picks lags, p = 0.  The unit root test produces a t-value of -1.51 which is greater than our 5% criterion -3.54.  Thus, at this point, we cannot reject the null of a unit root.

'Now, let's check for white noise.  To do that, I first set all the residuals = 0, then run the adf test and finally will check for white noise.

genr resid = 0
freeze(mode=overwrite,gexp_adf) gexp.uroot(adf,const,trend,info=sic)
freeze(gexp_adf_correl) resid.correl
 
'Based on the Q-statistic, I conclude that the residuals are white noise.  Putting it all together, I conclude that the gexp series is not level stationary.

'The next thing I do is test whether the differenced series is stationary using the ADF test.  I once again begin by graphing the (differenced) series.
 
genr gexpdiff = d(gexp)
freeze(figure_gexpdiff) gexpdiff.line
figure_gexpdiff.addtext(t) Dgexp (Pakistan):  1973-2008
figure_gexpdiff.addtext(b) Year
figure_gexpdiff.addtext(l) Dgexp
figure_gexpdiff.legend(off)

'The graph is not particularly illuminating.  Depending on how you look at it, it could have a time trend to it.  However, since the previous ADF test concluded that the equation for gexp do have a time trend, we run the ADF test with both a constant as well as a time trend.

'So we begin the whole process over again: 

genr gexpdiff = d(gexp)
freeze(table_9_8_gexpdiff1_adf) gexpdiff.uroot(adf,const,info=sic)

'Note that the SIC gexptomatic lag selection picks no lags, p =0.  The unit root test produces a t-value of -5.47 which is now smaller than our 5% criterion -2.95.  Thus, we may now reject the null of nonstationarity in first differenced series of gexp.  There is no reason to go further.  The last thing we do is check for white noise.

genr resid = 0
freeze(mode=overwrite,gexpdiff1_adf) gexpdiff.uroot(adf,const,info=sic)
freeze(gexpdiff1_adf_correl) resid.correl

''Based on the Q-statistic, I conclude that the residuals are white noise.  Putting it all together, I conclude that the gexp series is I(1).

'*********************************************************************************
'DF-GLS Unit Root Test for Pakistan's gexp
'We now run our first dfgls test
'**********************************************************************************
 
freeze(table_9_8_gexp_dfgls) gexp.uroot(dfgls,trend,info=sic)

'Note that the SIC automatic lag selection picks lags, p = 0.  The unit root test produces a t-value of -1.55 which is smaller than our 5% criterion -3.19.  Thus, at this point, we may reject the null of a unit root.

'Now let's see if the series is difference stationary or not

genr gexpdiff = d(gexp)
freeze(table_9_8_gexpdiff1_dfgls) gexpdiff.uroot(dfgls,const,info=sic)

'Note that the SIC automatic lag selection picks no lags, p = 0.  The unit root test produces a t-value of -4.09 which is now smaller than our 5% criterion -1.95.  Thus, we may not reject the null of nonstationarity in first differenced series of gexp.  
 
''Putting it all together, I conclude that the gexp series is I(1), a finding compatable with my ADF test results.

'***************************************************************
'Graph for Pakistan's  gdppc
'***************************************************************
                                        
genr gdppc = gdppc
freeze(figure_gdppc) gdppc.line
figure_gdppc.addtext(t) gdppc (Pakistan):  1973-2008
figure_gdppc.addtext(b) Year
figure_gdppc.addtext(l) gdppc
figure_gdppc.legend(off)
                                                 
'We see from the FIGURE that gdppc has a time trend to it.  So we would include both an intercept and a time trend in our unit root regression equtions. 

'*********************************************************************************
'ADF Unit Root Test for Pakistan's gdppc
'We now run our first ADF test
'**********************************************************************************
 
freeze(table_9_8_gdppc_adf) gdppc.uroot(adf,trend,info=sic)

'Note that the SIC automatic lag selection picks lags, p = 1.  The unit root test produces a t-value of -1.64 which is greater than our 5% criterion -3.55.  Thus, at this point, we cannot reject the null of a unit root.

'Now, let's check for white noise.  To do that, I first set all the residuals = 0, then run the adf test and finally will check for white noise.

genr resid = 0
freeze(mode=overwrite,gdppc_adf) gdppc.uroot(adf,const,trend,info=sic)
freeze(gdppc_adf_correl) resid.correl
 
'Based on the Q-statistic, I conclude that the residuals are white noise.  Putting it all together, I conclude that the gdppc series is not level stationary.

'The next thing I do is test whether the differenced series is stationary using the ADF test.  I once again begin by graphing the (differenced) series.
 
genr gdppcdiff = d(gdppc)
freeze(figure_gdppcdiff) gdppcdiff.line
figure_gdppcdiff.addtext(t) Dgdppc (Pakistan):  1973-2008
figure_gdppcdiff.addtext(b) Year
figure_gdppcdiff.addtext(l) Dgdppc
figure_gdppcdiff.legend(off)

'The graph is not particularly illuminating.  Depending on how you look at it, it could have a time trend to it.  However, since the previous ADF test concluded that the equation for gdppc do have a time trend, we run the ADF test with both a constant as well as a time trend.

'So we begin the whole process over again: 

genr gdppcdiff = d(gdppc)
freeze(table_9_8_gdppcdiff1_adf) gdppcdiff.uroot(adf,const,info=sic)

'Note that the SIC gdppctomatic lag selection picks no lags, p =0.  The unit root test produces a t-value of -4.57 which is now smaller than our 5% criterion -2.95.  Thus, we may now reject the null of nonstationarity in first differenced series of gdppc.  There is no reason to go further.  The last thing we do is check for white noise.

genr resid = 0
freeze(mode=overwrite,gdppcdiff1_adf) gdppcdiff.uroot(adf,const,info=sic)
freeze(gdppcdiff1_adf_correl) resid.correl

''Based on the Q-statistic, I conclude that the residuals are white noise.  Putting it all together, I conclude that the gdppc series is I(1).

'*********************************************************************************
'DF-GLS Unit Root Test for Pakistan's gdppc
'We now run our first dfgls test
'**********************************************************************************
 
freeze(table_9_8_gdppc_dfgls) gdppc.uroot(dfgls,trend,info=sic)

'Note that the SIC automatic lag selection picks lags, p = 1.  The unit root test produces a t-value of -1.63 which is smaller than our 5% criterion -3.19.  Thus, at this point, we may reject the null of a unit root.

'Now let's see if the series is difference stationary or not

genr gdppcdiff = d(gdppc)
freeze(table_9_8_gdppcdiff1_dfgls) gdppcdiff.uroot(dfgls,const,info=sic)

'Note that the SIC automatic lag selection picks no lags, p = 0.  The unit root test produces a t-value of -4.09 which is now smaller than our 5% criterion -1.95.  Thus, we may not reject the null of nonstationarity in first differenced series of gdppc.  
 
''Putting it all together, I conclude that the gdppc series is I(1), a finding compatable with my ADF test results.

'*********************************************
'Single Equation Cointegration Methods
'*********************************************

'******************************************************************
''Graph the suspected cointegrated series together
'******************************************************************
'The first step is to print out a graph of the series.  This is very important!

group g1 rer_def_NT rer_def_T rop gdppc gexp a_tilde
freeze(figure1) g1.line(x)
figure1.setelem(1) lcolor(black) symbol(1) lpat(1)
figure1.setelem(2) lcolor(black) symbol(4) lpat(1)
figure1.setelem(3) lcolor(black) symbol(7) lpat(1)
figure1.setelem(3) lcolor(black)
figure1.options linepat
figure1.addtext(t) rer_def_NT, lnpT & a_tilde (Pakistan & U.S): 1973-2008
figure1.addtext(b) Year
figure1.addtext(l) rer_def_NT
figure1.addtext(l) lnpT
figure1.addtext(r) a_tilde

''*******************************************************
''S1.A.Engle-Granger Approach to Cointegration
'*******************************************************

genr resid = 0
equation eg.ls rer_def_NT c rer_def_T rop gdppc gexp a_tilde
genr EC1 = resid

'First we test if the residuals of above regression are level stationary or not. If yes, next we'll proceed towards estimation of error correction model.

'***************************************************************
'Graph for Pakistan's  EC
'***************************************************************                                               
genr EC1 = EC1
freeze(figure_EC1) EC1.line
figure_EC1.addtext(t) EC1 (Pakistan):  1973-2008
figure_EC1.addtext(b) Year
figure_EC1.addtext(l) EC1
figure_EC1.legend(off)
                                                 
'We see from the FIGURE that EC has time trend to it.  So we would include both an intercept and trend in our unit root regression equtions. 

'*********************************************************************************
'EG Test for Cointegration
'**********************************************************************************
 
freeze(table_9_8_EGC) g1.coint(method=eg)

'The null hypothesis will not be rejected as suggested by sample statistics.

''******************************************
''S1.B.Error Correction Model (ECM)
'*******************************************
''*******************************************************
'Selecting the number of lags in the VAR  *
'*******************************************************
'NOTE: We do this because we need to have the "right" number of lags when it comes time to estimate our VEC model and test for cointegration.

var var1.ls 1 4   g1
freeze(var1_lagtest1) var1.laglen(4)
freeze(var1_lagtest2) var1.testlags

'The laglength test above indicates that the should have VAR has 4 lags. But we cannot afford to have this number of lags. 

var var2.ls 1 1  g1
freeze(var2_artest1) var2.correl
freeze(var2_artest2) var2.qstats(12)
freeze(var2_artest3) var2.arlm(12)

'The residuals are white noise.

'We now try different lags of d(rer_def_T) and d(a_tilde), comparing SIC values across specifications.

genr resid = 0
equation eg.ls rer_def_NT c rer_def_T rop gdppc gexp a_tilde
genr ec1 = resid

var table_9_8_eg2a.ls 0 0 d(rer_def_NT)   @  c ec1(-1) d(rer_def_NT(-1)) 

var table_9_8_eg2b.ls 0 0 d(rer_def_NT)   @  c ec1(-1)  d(rer_def_NT(-1)) d(rer_def_T(-1)) d(rop(-1)) d(gdppc(-1)) d(gexp(-1)) d(a_tilde(-1))

'The evidence suggests that Model A is best.  Now we test that model for serial correlation.

var table_9_8_eg2a.ls 0 0 d(rer_def_NT)   @   c ec1(-1) d(rer_def_NT(-1)) 
freeze(table_9_8_eg2a_artest1) table_9_8_eg2a.correl
freeze(table_9_8_eg2a_artest2) table_9_8_eg2a.qstats(12)
freeze(table_9_8_eg2a_artest3) table_9_8_eg2a.arlm(12)

'The residuals are absolutely white noise.

''*************************
''Estimating EC Model  
'**************************

'We'll now take the above specified model and turn it into an ECM. We shall run NW-HAC least squares model for establishing error correction mechanism.

'We now estimate the corresponding ECM:

equation table_9_8_ecm.ls(n) d(rer_def_NT) c ec1(-1) d(rer_def_NT(-1)) 

'Note that the SR effect is significant as the EC coefficient is of value -0.36 is statistically significant at better than 5% significance level.

''**********************************************
''S2.A & S2.B: Obtaining LR Coefficients
'***********************************************
'Now, by employing FMOLS and DOLS cointegration regression estimators, finally we shall calculate our LR coefficient i.e. BS coefficient for Pakistan against U.S.

equation table_9_8_LReqn1_fmols.cointreg(method=fmols) rer_def_NT rer_def_T rop gexp gdppc a_tilde

equation table_9_8_LReqn2_dols.cointreg(method=dols, trend=constant, lag=1,lead=1 ) rer_def_NT rer_def_T rop gexp gdppc a_tilde

'The BS coefficient obtained through FMOLS and DOLS estimators are statistically isignificant and/or bearing incorrect signs. Thus, there is invalid evidence in support of BS effect existing for Pakistan.

''******************************************
'Multivariate Cointegration Approach
'******************************************

''*************************************************************
''Check if the VAR(2) model is dynamically stable
'*************************************************************
freeze(table_9_8_var2_varstable) var2.arroots(graph)

'The model is dynamically stable.

''**********************************************************************
''M1.A & M1.B: Identifying the number of cointegrating vectors
'***********************************************************************
'Having identified the appropriate number of lags to put in, I now go on to test for the appropriate number of cointegrating equations.

freeze(table_9_8_var2_coint) var2.coint(s,1)

'This command estimates all possible combinations of constants and trends in the level data series and the cointegrating equations. All the results indicate 1 cointegrating vectors.
'
'GENERAL NOTE:, in practice, cases 1 and 5 are rarely used. One should use case 1 only if one knows that all series have zero mean. Case 5 may provide a good fit in-sample but will produce implausible forecasts out-of-sample. As a rough guide, use case 2 if none of the series appear to have a trend. For trending series, use case 3 if you believe all trends are stochastic; if you believe some of the series are trend stationary, use case 4.

'Note that the 5 cases are identified under "Johansen cointegration test" in Eviews. They run from most restrictive (no constants in either the level series or CEs) to most general (trend terms in both the level series and CEs).

''******************************************************************
''M2.A, M2.B & M3: Vector Error Correction Model (VECM)
'*******************************************************************

' For estimating the LR relationship, corresponding VEC command is:

var table_9_8_vec_c.ec(c,1)  0 0 rer_def_NT rer_def_T rop gexp gdppc a_tilde
var table_9_8_vec_d.ec(d,1)  0 0 rer_def_NT rer_def_T rop gexp gdppc a_tilde

'CONCLUSION:  I conclude that rer_def_NT and a_tilde are not cointegrated in the Pakistan's data.


